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Abstract By considering the fluctuation of grand potential Q around equilibrium with respect to small one-particle 
density fluctuations 6pq(7), the phase instability of restricted primitive model (RPM) of ionic systems is investigated. 
We use the integral equation theory to calculate the direct correlation functions in the reference hypernetted chain 
approximation and obtain the spinodal line of RPM. Our analysis explicitly indicates that the gas-fluid phase instability 
is induced by k = 0 fluctuation mode, while the fluid-solid phase instability is related to k 4 0 fluctuation modes. The 
spinodal line is qualitatively consistent with the result of computer simulations by others. 
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1 Introduction 

For ionic system the most widely investigated model is 
the restricted primitive model (RPM), in which the ions 
are modeled as equal-sized hard spheres of diameter a and 
carry charges of the same absolute value in a medium 
of dielectric constant D. The whole system is electri- 
cally neutral. A number of studies have been focused 
on the gas-fluid phase transition and the critical behav- 
ior of ionic liquids relying on RPM.!!! Comparatively less 
attention was paid to fluid-solid phase transition. Still- 
inger and Lovett!?] sketched the phase diagram includ- 
ing the fluid-solid phase transition for the charged hard 
sphere. In Ref. [3] Barrat used the density functional the- 
ory of freezing to calculate the fluid-solid phase diagram of 
RPM, and found that the solid changes from a disordered 
face-centered cubic (fcc) structure to a Caesium Chloride 
(CsCl) structure when the temperature decreases to be 
low enough. Monte Carlo (MC) computer simulation was 
applied to investigate the fluid-solid phase transition of 
RPM by the authors in Refs. [4,5]. Later Bresme et al. in 
Ref. [6] found a new ordered fcc structure at low temper- 
ature and high density by using MC simulation. The full 
coexistence curves of gas, fluid, and solid phases of RPM 
were found by means of computer simulation.|7] 

Since MC simulations mentioned above only gave the 
T-p phase diagram and did not present explicit analysis of 
how the control parameters (e.g. T, p) influence the phase 
transition characteristics. For this purpose, we focus on 
the gas-fluid and fluid-solid phase instability of RPM by 
means of soft mode analysis. Based on the correlation 
functions which are calculated in the framework of inte- 
gral equation theory with the reference hypernetted chain 
(RHNC) approximation, we present a complete discussion 
on the phase instability of the RPM and obtain the spin- 
odal line, which is the border line between stability and 
instability regime either of gas-fluid transition or of fluid- 
solid transition. Our results explicitly show that the gas- 


fluid phase instability is related to k = 0 fluctuation mode, 
while the fluid-solid phase instability is induced by k 4 0 
fluctuation modes. 


2 Model 

Progress has been made in the theories and simulations 
of the restricted primitive model (RPM) with equal num- 
ber of + ions interacting by coulombic force in a medium of 
dielectric constant D. All the ions are considered as hard 
spheres with the same diameter a and carrying charges 
da = +qo, qo is the charge of a proton. The ion-ion po- 
tential between ion a and ion ĝ is given by 


Oo, 
Vap(T) =) dade 


r<a; 


(1) 


3 Method for Analysis of Instability 

In order to study the stability of the equilibrium state 
for the binary mixture, we consider the fluctuation of the 
grand potential around the minimum caused by small one- 
particle density fluctuations 6p,(7) and expand the grand 
potential to the second order 


6Q = QUT, V, Has Pa + dpal7)| — Q(T, V, Mas Pa) 


1 672 
a I pa (Ti )dpa(F2) le. 


x Ôpa(T1)ópa (T2) FLAT, (2) 
where a, 8 = 1,2, and 1, 2 represent the two components 


of mixture. The second functional derivatives of Q are 
related to the direct correlation functions cag(1, 2)89] by 


62 6a80(71, T2) 
To. far. fae eee =— Ca I. 2 . 3 
AA) a OAR 
For cag(1,2), 1 = 7 and 2 = 7. For a homogeneous 
system, the Fourier transformation of Eq. (2) reads 
1 kpT 1 
l= >= > (k) = = 
2s (e) 2 V 
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ap — (Papa)’/*Gaa(k)]5p6(k) - (4) 


x SOS Salk) [6 
k aß 


Here the integral over k is replaced by a sum over discrete 
k-values: 


1 1 g Spa(k) 
ay | thoy dpalk) = FE ’ 
k a 


and alk) is the Fourier transformation of pa(r). It is 
more convenient to discuss the physical meaning of Eq. (4) 
by introducing two new fluctuation variables as follows 


bp(k) = p7"? [5px (k) + 6p2(k)] , (5) 
Selk) = p~*/? (cre) /? p281 (k) — prdp2(k)] (6) 
where p = pı + p2, C1 = pi/p = c, and c = p2/p = 1- c. 
The fluctuation variable df(k) corresponds to total one- 


particle density fluctuation and ôc(k) corresponds to the 
concentration fluctuation. Now 6 can be rewritten as 


bart g ‘ 
IOS y 2a 
Mop(k) Mock) \ ( 6p(k) 


where the matrix M(k) is defined as 


M,,(k) = 1 — pleféi(k) + ch@20(k) + 2c1c2č12(k)], (8) 
Mec(k) = 1 — peice|@11(k) + G22(k) — 2@12(k)] , (9) 
Myc(k) = Mep(k) = (c12) pleaéaa(k) (9) 

— c1@11(k) — (c2 — c1)€12(k))]. (10) 


The matrix M(k) can be diagonalized with eigenvalues 
Aa(k) and eigenvectors a(k), and then the grand poten- 
tial fluctuation (7) becomes a sum of pure squares 


kT ISSA CR)? + Aa(R)ADS(AI], (11) 


6Q = = 
2 ME 


where 
Opi (k) = xi1(k)ðp(k) + xi, 2(k)ðc(k) , 


teie T [Mpp(k) + Meol) 


F ((Mpp(k) — Mec(k))” + 4Mpc(k))"/7]. 2) 


The system is stable if the grand potential increases for 
any fluctuation of one-particle densities. The fluctuation 
of the grand potential takes a quadratic form in Eq. (11) 
and the system is stable only when the eigenvalues of the 
coefficient matrix are positive. As a result, the border of 
a stability region, the spinodal line, is reached when the 
smaller eigenvalue A;(k) goes to zero, and we can get the 
corresponding ôp; (k), which is most unstable. The eigen- 
vector £1(k) characterizes the phase transition uniquely. 
When Z1(k) is close to 7e = (0,1), the instability is dom- 
inated by demixing. If 7ı(k) is close to 7p = (1,0), the 
instability is dominated by condensation. 

For RPM the two components of the mixture 1, 2 are 
the cation and anion respectively. One has p} = p- = 
p/2,c, =c_ = 1/2, č}4 (k) = č- (k), and then Eqs. (8), 
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(9), and (10) become 
Myp(k) = 1- Zoë (k) +E), (18) 
Melk) =1—Solé++(k)—@4-(B)], 04) 
Mye(k) = Mep(k) =0. (15) 


Hence we can observe the pure condensation or demix- 
ing phase instability for RPM. When we change T and p, 
the pure condensation apy transition occurs if App(k) = 
Mop(k) = 0 and Acc(k) = Mec(k) £ 0 (i.e. App(k) goes to 
zero first), and the pure emia phase instability occurs 
if Aec(k) = 0 and App(k) Æ 0 (ie. Acc(k) goes to zero 
first). In order to calculate App(k) and Acc(k), we need to 
know the correlation functions č}4 (k) and é,_(k), which 
will be calculated in the next section. 


4 Calculation of Correlation Functions 

The correlation functions of ionic fluids can be cal- 
culated by integral equation theory which contains two 
equations.!!°l The first one is the Ornstein—Zernike (O.Z.) 
equation which gives the exact relation between the ion- 
ion total correlation functions hag(1,2) and the ion-ion 
direct correlation functions cyg(1,2). For homogeneous 
liquids, the O.Z. equation is 


haat 2) =Cap(Ls 2+ Py | hal, 8)cya(3,2)4%7, (16) 
gi 


where a, 8, y = + are the ionic species. This equation can 
be interpreted as that the total correlation of two ions 
is the sum of direct correlation and indirect correlation 
transferred by all other particles. 

For homogenous fluids the correlation functions de- 
pend on only the distance between two particles. We 
have then hag(1,2) = hag(r) and cag(1,2) = cag(r) with 
r = |r" = ril. 

Since the O.Z. equation contains two unknown func- 
tions hag(r) and cag(r), another relation between direct 
and total correlation functions is needed. This is 


1+ hoe(r) =exp[—Bvag(r)+ hag(r)— Caa(r) + bag(r)], (17) 
where vqg(r) is the interaction potential between a-ion 
and (-ion, bag(r) is the so-called bridge term which can 
not be obtained exactly. The bridge term can be expressed 


as 55 
= RYN) 
n=3 


where Ro (r) represents the contribution of all n- 
particle direct correlation functions. Different forms of 
bag(r) lead to different approximations. In the RHNC 
approximation,!!!!7] ba (r) is replaced by the bridge 
term of a reference system (usually the hard sphere 
system), |13-14] 


baa (r) = nfl +hag (r)] hag (r) +g (r), (19) 
where hl? (r) and cf/?(rr) are the total and direct correla- 
tion functions of hard spheres, which can be calculated by 


(18) 


r>a, 
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using Verlet and Weis’s method.!!3:!4] The MC simulation 
for the ion-dipole mixture has shown that the RHNC ap- 
proximation works very well!!®) for ionic fluids. So we take 
the RHNC approximation here to calculate the correlation 
functions. 

The integral equation including O.Z. equation (16) 
and RHNC equation (17) can be solved numerically. [0] 
Because of the long range coulombic interaction, the di- 
rect correlation function has also a long range tail, which 
should be treated analytically." The total correlation 
functions decay exponentially with the Debye screening 
length €, as the correlation length. At very low densities 
Ep can be very large and the total correlation functions de- 
cay then not so fast. According to the standard procedure 
to deal with long range interactions, in our calculation, we 
divide the correlation functions into two parts: the short 
range part and the long range part which are dealt with 
analytically.[6] The short range part can be calculated by 
numerical iteration method. 


5 Results for RPM 

Using the ion-ion direct correlation functions cag(r) 
obtained from the RHNC integral equation, we calculate 
Mop(k) and Mec(k) in Eqs. (13) and (14) for different tem- 
peratures and densities. Usually, only the homogeneous 
fluctuation mode k = 0 is considered. Figure 1 shows 
App(k = 0) for the reduced temperature T* at p* = 0.04. 
Here T* = kpTDa/q and p* = pa? denote the reduced 
temperature and the reduced density. Figure 1 shows that 
App(k = 0) decreases as the decrease of the reduced tem- 
perature, while Acc(k) diverges when k — 0. It means 
App(k = 0) approaches to zero first as the decrease of T* 
for the fixed reduced temperature p*. The reduced tem- 
perature corresponding to the gas-fluid phase instability 
point T is obtained through extrapolating the values of 
App(k = 0) to zero. For example, see the solid curve in 
Fig. 1. 


App (k = 0) 


p* = 0.04 


0.08 0.12 
Beg 


0.16 0.20 


0.04 


Fig. 1 App (k = 0) for the reduced temperature at the 
reduced density p* = 0.04 in RHNC approximation. The 
solid curve is the result of extrapolation, and TŽ is the 
reduced temperature on spinodal line. 
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If App(k) and Acc(k) for k # 0 are considered, we can 
also study fluid-solid phase instability. In this article we 
study it. The curves in Fig. 2 respectively correspond to 
App(k) and àcc(k) in RHNC approximation. 

Figure 2 shows that the minimum values of App(k) are 
always smaller than that of Acc(k) for different reduced 
densities, and the minimum value of A,,(k) approaches 
to zero first as the increase of the reduced density for 
the fixed reduced temperature T* = 1.0. To make it 
clearer, we summarize the minimum values of App and 
Acc for different reduced densities p* in Fig. 3. The re- 
duced density corresponding to the instability point p% 
is obtained through extrapolating the minimum values 
of App(k) to zero. A condensation phase transition hap- 
pens at the value of reduced density pš for T* = 1.0. 
We also notice that App(k) always goes to zero first. It’s 
surely not a demixing transition (for which ,.(k) goes 
to zero first) since the same type of ions cannot assem- 
ble together. In addition, Fig. 2 shows that the minimum 
value of App(k) occurs at k # 0. Unlike the homogeneous 
phase instability indicated by k = 0, it denotes the fluid- 
solid phase instability caused by the density fluctuations. 
pe = 1.36 < V2 = px,,x, and p*,,, is the close packing 
reduced density, which implies our result is reasonable. 


2.0 


1.6 
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0.8 


0.4 


0.0 


Fig. 2 The comparison between àpp(k) and Acc(k) in 
RHNC approximation. T* = 1.0, p* from 0.4 to 1.1 with 
increment of 0.1. 


For T* = 0.05, which is the reduced temperature 
lower than gas-fluid critical reduced temperature, we show 
App(k) for different reduced densities in Fig. 4. The mini- 
mum values of \,,(k) occur at k 4 0 for p* = 0.4, 0.5, 0.6; 
while at k = 0 for p* = 0.22, 0.26. 

The minimum values of A,,(k) for different reduced 
density at T* = 0.05 is shown in Fig. 5. For p* < 0.3 the 
minimum values of App(k) occur at k = 0. The minimum 
value of App(k) decreases as the reduced density decreases, 
and when it goes to zero, a homogeneous phase instability, 
i.e. gas-fluid phase instability happens. For p* > 0.3, the 
minimum values of \,,(k) occur at k # 0, and with the 


No. 2 


reduced density increasing the minimum value of App(k) 
decreases and approaches to zero. At the point of A =0 
an inhomogeneous phase instability, i.e. fluid-solid phase 
instability, occurs. 


0.8 m m i Amin T*=1.0 
cc 


Fig. 3 The minimum value of A,,(k) and Acc(k) as the 
functions of the reduced density p* in RHNC approxi- 
mation at the reduced temperature T* = 1.0. The solid 
curve is the result of extrapolation, and p% is the reduced 
density on spinodal line. 
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Fig. 4 A,,(k) at reduced temperature T* = 0.05 and 
reduced densities p* = 0.22, 0.26, 0.4, 0.5, 0.6 respectively 
in RHNC approximation. 


Using the methods mentioned above we get the spin- 
odal line of RPM. See Fig. 6. The pš < 0.4 region corre- 
sponds to gas-fluid phase instability and p% > 1.0 region 
corresponds to fluid-solid phase instability. The RHNC 
integral function theory cannot give the result of the low 
temperature region between gas-fluid spinodal line and 
fluid-solid spinodal line possibly due to the complexity of 
triphasic coexistence. Comparing our result with the com- 
puter simulation results (Fig. 6 in Ref. [7]), we find that 
our result of gas-fluid phase transition fits the simulation 
result very well, but for fluid-solid phase transition our 
result is only qualitatively consistent with the simulation. 
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It reveals that the RHNC approximation does not work 
well for high reduced densities. 
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Fig. 5 The minimum values of A,,(k) as the function 
of the reduced density p* with the reduced temperature 
T* =0.05. 
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Fig. 6 The spinodal curve for the RPM in RHNC ap- 
proximation. 


6 Conclusions 

In this paper we use the integral equation theory in the 
reference hypernetted chain approximation to analyse the 
gas-fluid and fluid-solid phase instability of the restricted 
primitive model of ionic liquids. 

Our analysis is based on the correlation functions of 
the fluid mixture. The correlation functions are calcu- 
lated in the framework of integral equation theory. The 
integral equation contain two equations: the Ornstein- 
Zernike equation and the closure. RHNC approximation 
is applied as the closure approximation. We use the itera- 
tion method to solve the integral equation and then we can 
obtain the direct and total ion-ion correlation functions. 

We discuss the phase stability of a system by the 
grand potential fluctuation around equilibrium with re- 
spect to small one-particle density fluctuations in detail. 
The grand potential fluctuation is a quadratic form of 
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density fluctuations and the coefficient matrix can be ex- 
pressed by the direct correlation functions of the system. 
This matrix can be diagonalized by its eigenvalues and 
eigenvectors. Once one of the eigenvalues goes to zero, 
the system becomes unstable and the related eigenvector 
characterizes the phase instability uniquely. 

For RPM there are two kinds of eigenvalues A,,(k) and 
Acc(k). Using the correlation functions obtained from the 
RHNC integral equation theory, we can calculate App(k) 
and Acc(k) at different reduced densities and tempera- 
tures. If App(k) goes to zero first, it denotes a condensa- 
tion phase instability caused by the density fluctuations. 
If Acc(k) goes to zero first, it denotes a demixing phase 
instability caused by the concentration fluctuations. 

For a given reduced temperature and density, the min- 
imum value of App(k) is smaller than that of Acc(k), and 
the minimum value of App(k) goes to zero first. It means 
that the phase instability of RPM is a condensation phase 
transition caused by the density fluctuation. 
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Further investigation shows that if the minimum value 
of App(k) becomes zero at the fluctuation mode of k Æ 0, 
this phase instability is fluid-solid phase instability. On 
the other hand, if the minimum value of A,,(k) becomes 
zero at the fluctuation mode k = 0, this phase instabil- 
ity is a homogeneous condensation phase instability, i.e. 
gas-fluid phase instability. 

Under the integral equation theory the minimum value 
of App(k) cannot arrive at zero but very close. The spin- 
odal line of RPM is obtained by extrapolating the mini- 
mum value of A,,(k) to zero. Comparing our result with 
the phase diagram derived from computer simulation, we 
find that our result is qualitatively consistent with the re- 
sult of simulation, and at the low reduced densities where 
gas-fluid phase instability occurs, our result fits the simu- 
lation very well. 
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